home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 19 / CU Amiga Magazine's Super CD-ROM 19 (1998)(EMAP Images)(GB)[!][issue 1998-02].iso / CUCD / Programming / LEDA / source / src / plane / _simplex.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-16  |  13.3 KB  |  563 lines
  1. /*******************************************************************************
  2. +
  3. +  LEDA  3.1c
  4. +
  5. +
  6. +  _simplex.c
  7. +
  8. +
  9. +  Copyright (c) 1994  by  Max-Planck-Institut fuer Informatik
  10. +  Im Stadtwald, 6600 Saarbruecken, FRG     
  11. +  All rights reserved.
  12. + 
  13. *******************************************************************************/
  14.  
  15.  
  16. /****************************************************************************
  17. *                                                                           *
  18. *    simplex.c                                                              *
  19. *    *********                                                              *
  20. *    ist ein Programm zur Loesung linearer Programme !                      *
  21. *    *************************************************                      *
  22. *                                                                           *
  23. *    Die Implementierung basiert auf dem Zweiphasen-Simplex-Verfahren       *
  24. *    aus Dinkelbach, Operations Research, Kapitel 1.                        *
  25. *                                                                           *
  26. *    Kommentare und die Wahl der Bezeichnernamen beziehen sich in der       *
  27. *    Regel auf die entsprechenden Stellen im Buch.                          *
  28. *                                                                           *
  29. *    29.06.1992                                                             *
  30. *                                                                           *
  31. *    Sascha Portz und Johannes Helders                                      *
  32. *                                                                           *
  33. ****************************************************************************/
  34.  
  35.  
  36. #include <LEDA/simplex.h>
  37.  
  38.  
  39. static void Init_Tableau (matrix A, matrix& M, vector c, vector b,int kl,int gr,
  40.            int gl, int *&B, int *&NB, int zeilen, int spalten, int n)
  41. {
  42.    int nb = n;                      // Zaehler zur Initialisierung der Menge NB
  43.  
  44.    for ( int s=0 ; s<n; s++ )
  45.       M(0,s+1) = -c[s];
  46.  
  47.    for (int z=0; z<kl; z++ )
  48.    {
  49.       for ( int s=0; s < n; s++ )
  50.      M(z+1,s+1) = A(z,s);
  51.       M(z+1,n+z+1) = 1;              // Schlupfvariablen setzen
  52.       B[z] = n+z+1;                  // der z.te Basisvektor
  53.    }
  54.  
  55.    for ( z=0; z<gr; z++ )
  56.    {
  57.       for ( int s=0; s < n; s++ )
  58.      M(kl+z+1,s+1) = A(kl+z,s);
  59.       M(z+kl+1,n+z+kl+1) = -1;        // bei >= Nebenbedingungen neg. Schlupf
  60.       M(z+kl+1,n+kl+gr+z+1) = 1;      // kuenstliche Variable
  61.       B[kl+z] = n+kl+gr+z+1;          // der kl+gr+z.te Basisvektor
  62.       NB[nb] = n+kl+z+1; nb++;        // n+kl+gr+z ist Nichtbasisvaiable
  63.       M(zeilen-1,n+kl+gr+z+1) = 1;    // Initialisieren der Hilfszielfunktion
  64.    }
  65.  
  66.    for ( z=0; z<gl; z++ )
  67.    {
  68.       for ( int s=0; s<n; s++ )
  69.      M(z+kl+gr+1,s+1) = A(z+kl+gr,s);
  70.       M(z+kl+gr+1,n+kl+2*gr+z+1) = 1;  // kuenstliche Variable
  71.       B[kl+gr+z] = kl+2*gr+n+z+1;      // der z.te Basisvektor
  72.       M(zeilen-1,n+kl+2*gr+z+1) = 1;   //Initialisieren der Hilfszielfunktion
  73.    }
  74.  
  75.    for ( z=0; z<b.dim(); z++)
  76.       M(z+1,0) = b[z];
  77.  
  78.    for( s=0; s < n ; s++)               // 1,...,n sind Nichtbasisvariablen
  79.       NB[s] = s+1;
  80.  
  81.    for ( z=kl+1; z<=kl+gr+gl; z++)      // Subtrahieren der Zeilen mit kuenstl.
  82.       for ( s=0; s<spalten; s++)       // Variablen von der letzten Zeile
  83.      M(zeilen-1,s) = M(zeilen-1,s) - M(z,s);
  84.  
  85. }
  86.  
  87.  
  88. static void InitIntArray (int *&A, int *B, int l)
  89.  
  90. // Initialisiert das int-Array A der Laenge l mit dem int-Array B
  91.  
  92. {
  93.    for (int i=0; i<l; i++)
  94.       A[i] = B[i];
  95. }
  96.  
  97.  
  98. static void InitMatrix (matrix &A, matrix B, int dim1, int dim2)
  99.  
  100. // Initialisierung der Matrix A mit der Matrix B
  101.  
  102. {
  103.    for (int l=0; l<dim1; l++)
  104.       for (int m=0; m<dim2; m++)
  105.      A(l,m) = B(l,m);
  106. }
  107.  
  108.  
  109. static int PivotSpalte (matrix M, int *&NB, int z, int l, int m)
  110.  
  111. // in der Zeile z (z=0 oder z = m+1) wird die Pivotspalte gesucht
  112.  
  113. // l ist die Groesse der Menge der Nichtbasisvariablen
  114.  
  115. // accept = 1 impliziert die Gueltigkeit der gefundenen Pivotspalte, d.h.
  116. // es gibt in dieser Zeile ein Element, das > 0 ist,
  117. // ist dies nicht der Fall so wird nach einer neuen Pivotspalte gesucht
  118. // iter zaehlt die Anzahl der Iterationen
  119.  
  120. {
  121.    double eps = 0.00001;    // das Intervall [-eps..eps] entspricht NULL
  122.  
  123.    double min=-eps;
  124.    double minmin=-10000;
  125.    int mins = 0;
  126.    int p, accept = 0;
  127.    int iter = 1;
  128.  
  129.    while ( (accept == 0) && (iter <= l) )
  130.    {
  131.       for (int q=0; q<l; q++)
  132.      if ( ( M(z,NB[q]) <= min ) && ( M(z,NB[q]) > minmin ) )
  133.      {
  134.         min = M(z,NB[q]);
  135.         mins = NB[q];
  136.      }
  137.  
  138.       p = 1;                                // Test ob die gefundene Spalte
  139.       while ( (p <= m) && (accept == 0) )   // als Pivotspalte akzeptiert
  140.       {                                     // wird
  141.      if ( M(p,mins) > eps )
  142.         accept = 1;
  143.      p++;
  144.       }
  145.  
  146.       if (accept == 0)                      // falls nicht
  147.       {                                     // -> Suche nach der naechst-
  148.      minmin = min;                     // besten Pivotspalte
  149.      min = -eps;
  150.      q = 0;
  151.      mins = 0;
  152.       }
  153.       iter++;
  154.    }
  155.  
  156.    return mins;
  157. }
  158.  
  159.  
  160. static int D_PivotSpalte (matrix M, int *&NB, int z,int n,int kl,int gr)
  161.  
  162. // Liefert die Pivotspalte fuer einen Dualsimplexschritt in der Zeile z
  163.  
  164. {
  165.    double eps = 0.00001;    // das Intervall [-eps..eps] entspricht NULL
  166.  
  167.    double min = 1000;
  168.    int mins = 0;
  169.  
  170.    for(int q=0; q<n+gr; q++)
  171.       if ( ( M(z,NB[q]) < -eps ) && ( NB[q] <= n+kl+gr ) )
  172.      if ( min >= -M(0,NB[q]) / M(z,NB[q]) )
  173.      {
  174.         min = -M(0,NB[q]) / M(z,NB[q]);
  175.         mins = NB[q];
  176.      }
  177.    return mins;
  178. }
  179.  
  180.  
  181. static int PivotZeile (matrix M, int s, int m)
  182.  
  183. // Liefert die Pivotzeile fuer einen Primalsimplexschritt in der Spalte s
  184.  
  185. {
  186.    double eps = 0.00001;    // das Intervall [-eps..eps] entspricht NULL
  187.  
  188.    int minz = 0;
  189.    double min = 1000;
  190.  
  191.    for (int p=1; p<=m; p++)
  192.       if ( M(p,s) > eps )
  193.       if (min > M(p,0)/M(p,s) )
  194.       {
  195.          min = M(p,0)/M(p,s);
  196.          minz = p;
  197.       }
  198.    return minz;
  199. }
  200.  
  201.  
  202. static void Basiswechsel (matrix &M, int p, int q, int n,int kl, int gr, int gl)
  203.  
  204. // Basiswechsel : siehe Dinkelbach
  205.  
  206. {
  207.    double pivot = M(p,q);
  208.    double pivot2;
  209.    int spalten = n+kl+2*gr+gl;
  210.    int zeilen = kl+gr+gl;
  211.  
  212.    for (int s=0; s<=spalten; s++)
  213.       M(p,s)= M(p,s)/pivot;
  214.    for (int z=0; z<p; z++)
  215.    {
  216.       pivot2 = M(z,q);
  217.       for (int s=0; s<=spalten; s++)
  218.      M(z,s) = M(z,s)-M(p,s)*pivot2;
  219.    }
  220.    for (z=p+1; z<=zeilen+1; z++)
  221.    {
  222.       pivot2 = M(z,q);
  223.       for (int s=0; s<=spalten; s++)
  224.      M(z,s) = M(z,s)-M(p,s)*pivot2;
  225.    }
  226. }
  227.  
  228.  
  229. static void NB_NNB_Update (int *&B, int *&NB, int p, int q, int l)
  230.  
  231. // Nach jedem Simplex-Schritt werden die Mengen der Basis- ( B ) und
  232. // Nichtbasisvariablen ( NB ) neu berechnet
  233.  
  234. // p und q bezeichnen die Pivotzeile bzw die Pivotspalte
  235.  
  236. // l ist die Groesse der Menge der Nichtbasisvariablen (diese aendert
  237. // sich im Laufe des Verfahrens)
  238.  
  239. {
  240.    int h1, h2;
  241.  
  242.    for (int x=0; x<l; x++)
  243.       if (NB[x] == q)
  244.      h1 = x;
  245.    h2 = B[p-1];
  246.    B [p-1] = q;
  247.    NB [h1] = h2;
  248. }
  249.  
  250.  
  251. static void NNB_Update (int *&NB, int n, int kl, int gr, int gl)
  252.  
  253. // Veringerung der Menge der Nichtbasisvariablen um die k’nstlichen Variablen
  254.  
  255. {
  256.    int *nnb = new int[n-gl];
  257.    int m = 0;
  258.  
  259.    for (int l=0; l<n+gr; l++)
  260.       if ( NB[l] <= n+kl+gr )
  261.       {
  262.      nnb[m] = NB[l];
  263.      m++;
  264.       }
  265.    delete NB;
  266.    NB = nnb;
  267. }
  268.  
  269.  
  270. static int nicht_optimal (matrix &M, int n, int kl, int gr)
  271.  
  272. // returniert 1, wenn das Simplextableau nicht optimal ist, sonst 0
  273.  
  274. {
  275.    double eps = 0.00001;    // das Intervall [-eps..eps] entspricht NULL
  276.  
  277.    int not_opt = 0;
  278.    for (int i=1; i<=n+gr+kl; i++)
  279.       if ( M(0,i) < -eps )
  280.      not_opt = 1;
  281.    return not_opt;
  282. }
  283.  
  284.  
  285. static void Eliminiere_kuenstl_Var (matrix M, int *&B, int *&NB, int &error,
  286.                  int m, int n, int kl, int gr, int gl)
  287.  
  288. // die kuenstlichen Variablen, die noch Basisvariablen sind, werden zu
  289. // Nichtbasisvariablen "pivotiert"
  290.  
  291. {
  292.    int q;
  293.  
  294.    int l=0;
  295.    while ( (l<m) && nicht_optimal(M,n,kl,gr) )
  296.    {
  297.       if (B[l] > n+kl+gr)
  298.       {
  299.      q = D_PivotSpalte (M,NB,l+1,n,kl,gr);
  300.          if ( (q != 0) && ( error == 0 ) )
  301.      {
  302.         Basiswechsel (M,l+1,q,n,kl,gr,gl);
  303.         NB_NNB_Update (B,NB,l+1,q,n+gr);
  304.      }
  305.      else
  306.         error = 1;
  307.       }
  308.       l++;
  309.    }
  310.    if ( (error==0) && nicht_optimal(M,n,kl,gr) )
  311.       NNB_Update (NB,n,kl,gr,gl);
  312. }
  313.  
  314.  
  315. static vector Basisloesung (matrix M, int *B, int m, int n)
  316.  
  317. // returniert eine optimale Basisloesung
  318.  
  319. {
  320.    vector X(n);
  321.  
  322.    for (int l=0; l<m; l++)
  323.       if ( B[l] <= n )
  324.      X[B[l]-1] = M(l+1,0);
  325.  
  326.    return X;
  327. }
  328.  
  329.  
  330. static vector Extremalstrahl (matrix M, int *&B, int *&MF, int l, int m, int n)
  331.  
  332. // returniert die Steigung des Extremalstrahls
  333.  
  334. {
  335.    vector s(n);
  336.  
  337.    for (int j=0; j<m; j++)
  338.       if ( B[j] <= n )
  339.     s[B[j]-1] = -M(j+1,MF[l]);
  340.  
  341.    if ( MF[l] <= n )
  342.       s[MF[l]-1]=1;
  343.  
  344.    return s;
  345. }
  346.  
  347.  
  348. static int ZERO (vector v,int n)
  349.  
  350. // returniert 1, falls v gleich dem Nullvektor ist, sonst 0
  351.  
  352. {
  353.    double eps = 0.00001;    // das Intervall [-eps..eps] entspricht NULL
  354.  
  355.    int zero = 1;
  356.  
  357.    for (int i=0; i<n; i++)
  358.       if ( (v[i] < -eps) || (v[i] > eps) )
  359.          zero = 0;
  360.  
  361.    return zero;
  362. }
  363.  
  364.  
  365. static void Init_Matrix (matrix &A, vector x, int n)
  366.  
  367. // initialisiert die (1xn)-Matrix mit dem Vektor x
  368.  
  369. {
  370.    for (int i=0; i<n; i++)
  371.       A(0,i) = x[i];
  372. }
  373.  
  374.  
  375. static void Enlarge_Matrix (matrix &A, vector v, int n)
  376.  
  377. // vergroessert die Matrix A um den Zeilenvektor v
  378.  
  379. {
  380.    int dim1 = A.dim1();
  381.    matrix B(dim1+1,A.dim2());
  382.  
  383.    InitMatrix (B,A,dim1,A.dim2());
  384.  
  385.    for (int i=0; i<n; i++)
  386.       B(dim1,i) = v[i];
  387.  
  388.    A = B;
  389. }
  390.  
  391.  
  392. static void Test_auf_Extremalstrahl (matrix &M,int *&B,int *&NB,list<matrix>& L,
  393.                   int *&MF, int *&PZ, int i, vector x, int m, 
  394.                               int n, int kl, int gr)
  395.  
  396. // Berechnung und Ausgabe eines optimalen Extremalstrahls
  397. // Aufbau der Liste L der optimalen Loesungen
  398.  
  399. {
  400.    vector s(n);
  401.    matrix A (1,n);
  402.  
  403.    Init_Matrix (A,x,n);
  404.  
  405.    for (int l=0; l<i; l++)
  406.       if ( NB[l] <= n+gr+kl )             // Test ob im optimalen Tableau ein
  407.      if ( M(0,NB[l]) == 0 )        // Zielfunktionskoeffient einer Nicht-
  408.         MF[l] = NB[l];                // basisvariable = 0 ist
  409.      else
  410.         MF[l] = 0;
  411.       else
  412.          MF[l] = 0;
  413.  
  414.    for ( l=0; l<i; l++)
  415.    {
  416.       if ( MF[l] > 0 )
  417.       {
  418.      PZ[l] = PivotZeile (M,MF[l],m);
  419.      if ( PZ[l] == 0 )
  420.      {
  421.         s = Extremalstrahl (M,B,MF,l,m,n);
  422.         if ( !ZERO(s,n) )
  423.            Enlarge_Matrix (A,s,n);
  424.      }
  425.       }
  426.    }
  427.  
  428.    L.append (A);
  429. }
  430.  
  431.  
  432. static
  433. void Test_auf_Mehrfachloesungen(matrix &M,int *&B,int *&NB,list<matrix> &L,
  434.                 int x, vector v,int zeilen, int spalten,
  435.                 int m, int n, int kl, int gr, int gl)
  436.  
  437. // Berechnung von eventuell existierenden weiteren optimalen Basisloesungen
  438. // x ist die Groesse der Menge der Nichtbasisvariablen
  439.  
  440. {
  441.                     // Hilfsmatrix zur Berechnung weiterer
  442.    matrix A (m+2,n+kl+2*gr+gl+1);   // optimaler Simplextableaus
  443.  
  444.    int *MF = new int[x];              // fuer die Indizes der Spalten in
  445.                       // denen die 0.te Zeile =0 ist
  446.    int *PZ = new int[x];              // die zugehoerigen Pivotzeilen
  447.  
  448.    int *H1 = new int[x];              // Hilfsarrays
  449.    int *H2 = new int[x];
  450.  
  451.    int *b = new int [m];              // Hilfsarrays fuer B
  452.    int *nb = new int [x];             // und NB
  453.  
  454.  
  455.    Test_auf_Extremalstrahl (M,B,NB,L,MF,PZ,x,v,m,n,kl,gr);
  456.  
  457.    for (int l=0; l<x; l++)
  458.    {
  459.       if ( MF[l] > 0 )
  460.      if ( PZ[l] > 0 )
  461.      {
  462.         InitMatrix (A,M,zeilen,spalten);
  463.         InitIntArray (b,B,m);
  464.         InitIntArray (nb,NB,x);
  465.         Basiswechsel (A,PZ[l],MF[l],n,kl,gr,gl);
  466.         NB_NNB_Update (b,nb,PZ[l],MF[l],x);
  467.         v = Basisloesung (A,b,m,n);
  468.         Test_auf_Extremalstrahl (A,b,nb,L,H1,H2,x,v,m,n,kl,gr);
  469.                          // H1, H2 sind Hilfsarrays
  470.      }
  471.    }
  472.    delete MF; delete PZ; delete H1;
  473.    delete H2; delete b; delete nb;
  474. }
  475.  
  476.  
  477. static
  478. void Run (matrix &M, int *&B, int *&NB, list<matrix> &L, int zeilen, 
  479.           int spalten, int m, int n, int kl, int gr, int gl)
  480.  
  481. // Hauptprozedur
  482.  
  483. {
  484.    double eps = 0.00001;    // das Intervall [-eps..eps] entspricht NULL
  485.  
  486.    int p, q;
  487.    int error=0;
  488.    vector x;
  489.  
  490.    q = PivotSpalte (M,NB,m+1,n+gr,m);
  491.    while (q != 0)
  492.    {
  493.       p = PivotZeile (M,q,m);
  494.       if (p != 0)
  495.       {
  496.      Basiswechsel (M,p,q,n,kl,gr,gl);
  497.      NB_NNB_Update (B,NB,p,q,n+gr);
  498.       }
  499.       q = PivotSpalte (M,NB,m+1,n+gr,m);
  500.    }
  501.  
  502.    if ( ( M(m+1,0) > -eps ) && ( M(m+1,0) < eps ) )
  503.    {
  504.  
  505.       Eliminiere_kuenstl_Var (M,B,NB,error,m,n,kl,gr,gl);
  506.       if ( nicht_optimal (M,n,kl,gr) )
  507.       {
  508.  
  509. // falls error = 1 so konnte ein notwendiger Dualsimplexschritt nicht
  510. // durchgefuehrt werden
  511.  
  512.      if ( error == 0 )
  513.      {
  514.         q = PivotSpalte (M,NB,0,n-gl,m);
  515.         while (q != 0)
  516.         {
  517.            p = PivotZeile (M,q,m);
  518.            if (p != 0)
  519.            {
  520.           Basiswechsel (M,p,q,n,kl,gr,gl);
  521.           NB_NNB_Update (B,NB,p,q,n-gl);
  522.            }
  523.            q = PivotSpalte (M,NB,0,n-gl,m);
  524.         }
  525.         x = Basisloesung(M,B,m,n);
  526.         Test_auf_Mehrfachloesungen(M,B,NB,L,n-gl,x,zeilen,
  527.                        spalten,m,n,kl,gr,gl);
  528.      }
  529.       }
  530.       else                              // das Tableau war schon optimal
  531.       {
  532.       x = Basisloesung(M,B,m,n);
  533.       Test_auf_Mehrfachloesungen (M,B,NB,L,n+gr,x,zeilen,
  534.                       spalten,m,n,kl,gr,gl);
  535.       }
  536.    }
  537. }
  538.  
  539.  
  540. list<matrix> SIMPLEX(matrix A, int i, int j, int k, vector b, vector c)
  541. {
  542.    int *B; int *NB;
  543.  
  544.    list<matrix> L;
  545.  
  546.    int kl =i,  gr=j, gl=k;
  547.    int n = A.dim2(),m = A.dim1();        // setzen der globalen Variablen
  548.    int zeilen = m+2, spalten = n+m+gr+1;
  549.  
  550.    matrix M(zeilen,spalten);
  551.    B  = new int[m];
  552.    NB = new int[n+gr];
  553.  
  554.    Init_Tableau(A,M,c,b,i,j,k,B,NB,zeilen,spalten,n);
  555.  
  556.    Run (M,B,NB,L,zeilen,spalten,m,n,kl,gr,gl);
  557.  
  558.    delete B; delete NB;
  559.  
  560.    return L;
  561. }
  562.  
  563.